Import data

# Read in bpi data
bpi <- readr::read_rds('./data/bpi.rds') %>%
    arrange(ID) %>%
    mutate(ID = stringr::str_to_lower(ID))

# Read in demographics.rds for site, group, and sex info
foo <- readr::read_rds('./data/demographics.rds') %>%
    select(ID, Site, Group, Sex) %>%
    rename(ID2 = ID) %>%
    arrange(ID2)

# Join the two datasets 
bpi <- foo %>%
    bind_cols(bpi) %>%
    rowwise() %>%
    mutate(ID_match = grepl(pattern = ID2,
                            x = ID))

# Check for mismatches from join
nrow(filter(bpi, ID_match == FALSE))

# Clean-up columns
bpi <- bpi %>%
    select(-ID,
           -ID_match) %>%
    rename(ID = ID2)

Quick look

glimpse(bpi)
## Observations: 160
## Variables: 88
## $ ID             <chr> "j1", "j10", "j11", "j12", "j17", "j18", "j19",...
## $ Site           <chr> "u1", "u1", "u1", "u1", "u1", "u1", "u1", "u1",...
## $ Group          <chr> "p", "t", "p", "t", "t", "t", "t", "p", "t", "t...
## $ Sex            <chr> "female", "male", "female", "female", "female",...
## $ Pain_Pres_BL   <chr> "Yes", "Yes", "Yes", "Yes", "Yes", NA, "Yes", "...
## $ Pain_Pres_Wk4  <chr> NA, "Yes", "Yes", NA, NA, NA, "Yes", NA, "Yes",...
## $ Pain_Pres_Wk8  <chr> NA, "Yes", "Yes", "Yes", NA, NA, "Yes", "Yes", ...
## $ Pain_Pres_Wk12 <chr> NA, "Yes", "Yes", "Yes", "Yes", NA, NA, NA, "Ye...
## $ Pain_Pres_Wk24 <chr> NA, "Yes", NA, "Yes", NA, NA, NA, "Yes", "Yes",...
## $ Pain_Pres_Wk48 <chr> NA, "Yes", NA, "Yes", NA, NA, NA, "Yes", NA, "Y...
## $ Worst_BL       <dbl> 8, 10, 9, 0, 10, NA, 9, 8, 10, 6, NA, 9, 10, 6,...
## $ Worst_Wk4      <dbl> NA, 9, 8, NA, NA, NA, 8, NA, 10, 8, NA, NA, NA,...
## $ Worst_Wk8      <dbl> NA, 9, 10, 4, NA, NA, 9, 6, 7, 7, NA, NA, NA, N...
## $ Worst_Wk12     <dbl> NA, 9, 10, 7, 10, NA, NA, NA, 7, 8, NA, NA, NA,...
## $ Worst_Wk24     <dbl> NA, 9, NA, 6, NA, NA, NA, 9, 5, 8, NA, NA, NA, ...
## $ Worst_Wk48     <dbl> NA, 8, NA, 8, NA, NA, NA, 6, NA, 8, NA, NA, NA,...
## $ Least_BL       <dbl> 4, 5, 3, 0, 3, NA, 3, 3, 2, 3, NA, 1, 4, 2, 2, ...
## $ Least_Wk4      <dbl> NA, 4, 2, NA, NA, NA, 3, NA, 3, 3, NA, NA, NA, ...
## $ Least_Wk8      <dbl> NA, 4, 5, 1, NA, NA, 4, 3, 3, 3, NA, NA, NA, NA...
## $ Least_Wk12     <dbl> NA, 6, 5, 4, 5, NA, NA, NA, 3, 3, NA, NA, NA, N...
## $ Least_Wk24     <dbl> NA, 5, NA, 2, NA, NA, NA, 7, 2, 3, NA, NA, NA, ...
## $ Least_Wk48     <dbl> NA, 5, NA, 2, NA, NA, NA, 2, NA, 4, NA, NA, NA,...
## $ Avg_BL         <dbl> 4, 7, 6, 0, 6, NA, 6, 6, 5, 5, NA, 4, 6, 4, 3, ...
## $ Avg_Wk4        <dbl> NA, 6, 5, NA, NA, NA, 5, NA, 5, 5, NA, NA, NA, ...
## $ Avg_Wk8        <dbl> NA, 7, 8, 2, NA, NA, 6, 5, 4, 5, NA, NA, NA, NA...
## $ Avg_Wk12       <dbl> NA, 7, 7, 6, 8, NA, NA, NA, 5, 5, NA, NA, NA, N...
## $ Avg_Wk24       <dbl> NA, 7, NA, 4, NA, NA, NA, 8, 3, 6, NA, NA, NA, ...
## $ Avg_Wk48       <dbl> NA, 7, NA, 5, NA, NA, NA, 4, NA, 6, NA, NA, NA,...
## $ Now_BL         <dbl> 6, 0, 9, 0, 3, NA, 0, 8, 6, 3, NA, 4, 10, 6, 1,...
## $ Now_Wk4        <dbl> NA, 4, 6, NA, NA, NA, 0, NA, 0, 5, NA, NA, NA, ...
## $ Now_Wk8        <dbl> NA, 5, 10, 1, NA, NA, 0, 7, 4, 3, NA, NA, NA, N...
## $ Now_Wk12       <dbl> NA, 5, 10, 3, 8, NA, NA, NA, 5, 8, NA, NA, NA, ...
## $ Now_Wk24       <dbl> NA, 8, NA, 3, NA, NA, NA, 9, 2, 8, NA, NA, NA, ...
## $ Now_Wk48       <dbl> NA, 2, NA, 5, NA, NA, NA, 2, NA, 2, NA, NA, NA,...
## $ Rx_BL          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ Rx_Wk4         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ Rx_Wk8         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ Rx_Wk12        <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ Rx_Wk24        <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ Rx_Wk48        <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ Relief_BL      <dbl> 0, 0, 30, 70, NA, NA, 70, 0, 40, 50, NA, 0, 10,...
## $ Relief_Wk4     <dbl> NA, 0, 30, NA, NA, NA, 40, NA, 80, 20, NA, NA, ...
## $ Relief_Wk8     <dbl> NA, 0, NA, 70, NA, NA, 70, NA, 60, 50, NA, NA, ...
## $ Relief_Wk12    <dbl> NA, 0, 10, 50, 10, NA, NA, NA, 50, 20, NA, NA, ...
## $ Relief_Wk24    <dbl> NA, 0, NA, 30, NA, NA, NA, NA, 70, 10, NA, NA, ...
## $ Relief_Wk48    <dbl> NA, 40, NA, 40, NA, NA, NA, NA, NA, 30, NA, NA,...
## $ Activ_BL       <dbl> 6, 0, 10, 0, 7, NA, 7, 7, 6, 4, NA, 9, 10, 5, 5...
## $ Activ_Wk4      <dbl> NA, 9, 0, NA, NA, NA, 2, NA, 2, 3, NA, NA, NA, ...
## $ Activ_Wk8      <dbl> NA, 8, 10, 4, NA, NA, 0, 6, 5, 0, NA, NA, NA, N...
## $ Activ_Wk12     <dbl> NA, 8, 10, 7, 10, NA, NA, NA, 5, 7, NA, NA, NA,...
## $ Activ_Wk24     <dbl> NA, 5, NA, 7, NA, NA, NA, 8, 5, 7, NA, NA, NA, ...
## $ Activ_Wk48     <dbl> NA, 6, NA, 7, NA, NA, NA, 7, NA, 8, NA, NA, NA,...
## $ Mood_BL        <dbl> 8, 8, 0, 0, 4, NA, 9, 10, 5, 0, NA, 10, 5, 10, ...
## $ Mood_Wk4       <dbl> NA, 6, 0, NA, NA, NA, 0, NA, 3, 5, NA, NA, NA, ...
## $ Mood_Wk8       <dbl> NA, 7, 3, 8, NA, NA, 0, 7, 3, 0, NA, NA, NA, NA...
## $ Mood_Wk12      <dbl> NA, 4, 10, 8, 10, NA, NA, NA, 4, 8, NA, NA, NA,...
## $ Mood_Wk24      <dbl> NA, 6, NA, 8, NA, NA, NA, 10, 5, 3, NA, NA, NA,...
## $ Mood_Wk48      <dbl> NA, 5, NA, 8, NA, NA, NA, 6, NA, 7, NA, NA, NA,...
## $ Walking_BL     <dbl> 4, 0, 7, 0, 8, NA, 0, 4, 6, 3, NA, 10, 9, 2, 0,...
## $ Walking_Wk4    <dbl> NA, 0, 0, NA, NA, NA, 4, NA, 4, 4, NA, NA, NA, ...
## $ Walking_Wk8    <dbl> NA, 7, 7, 0, NA, NA, 0, 7, 3, 4, NA, NA, NA, NA...
## $ Walking_Wk12   <dbl> NA, 7, 10, 0, 10, NA, NA, NA, 2, 10, NA, NA, NA...
## $ Walking_Wk24   <dbl> NA, 4, NA, 0, NA, NA, NA, 7, 5, 8, NA, NA, NA, ...
## $ Walking_Wk48   <dbl> NA, 4, NA, 3, NA, NA, NA, 3, NA, 8, NA, NA, NA,...
## $ Work_BL        <dbl> 5, 0, 8, 0, 6, NA, 5, 7, 6, 2, NA, 8, 10, 10, 6...
## $ Work_Wk4       <dbl> NA, 3, 0, NA, NA, NA, 4, NA, 4, 4, NA, NA, NA, ...
## $ Work_Wk8       <dbl> NA, 3, 10, 6, NA, NA, 0, 10, 3, 0, NA, NA, NA, ...
## $ Work_Wk12      <dbl> NA, 4, 10, 7, 8, NA, NA, NA, 4, 6, NA, NA, NA, ...
## $ Work_Wk24      <dbl> NA, 5, NA, 4, NA, NA, NA, 8, 1, 7, NA, NA, NA, ...
## $ Work_Wk48      <dbl> NA, 5, NA, 5, NA, NA, NA, 5, NA, 8, NA, NA, NA,...
## $ Rel_BL         <dbl> 4, 5, 0, 0, 7, NA, 9, 4, 5, 0, NA, 10, 6, 6, NA...
## $ Rel_Wk4        <dbl> NA, 0, 0, NA, NA, NA, 0, NA, 2, 0, NA, NA, NA, ...
## $ Rel_Wk8        <dbl> NA, 4, 6, 6, NA, NA, 0, 9, 3, 0, NA, NA, NA, NA...
## $ Rel_Wk12       <dbl> NA, 4, 7, 8, 7, NA, NA, NA, 2, 4, NA, NA, NA, N...
## $ Rel_Wk24       <dbl> NA, 6, NA, 5, NA, NA, NA, 10, 5, 3, NA, NA, NA,...
## $ Rel_Wk48       <dbl> NA, 6, NA, 6, NA, NA, NA, 5, NA, 7, NA, NA, NA,...
## $ Sleep_BL       <dbl> 2, 8, 10, 0, 0, NA, 0, 6, 0, 0, NA, 10, 10, 7, ...
## $ Sleep_Wk4      <dbl> NA, 6, 10, NA, NA, NA, 0, NA, 4, 0, NA, NA, NA,...
## $ Sleep_Wk8      <dbl> NA, 0, 10, 3, NA, NA, 6, 0, 3, 0, NA, NA, NA, N...
## $ Sleep_Wk12     <dbl> NA, 6, 10, 0, 10, NA, NA, NA, 0, 0, NA, NA, NA,...
## $ Sleep_Wk24     <dbl> NA, 4, NA, 0, NA, NA, NA, 10, 1, 1, NA, NA, NA,...
## $ Sleep_Wk48     <dbl> NA, 7, NA, 5, NA, NA, NA, 2, NA, 7, NA, NA, NA,...
## $ Enjoy_BL       <dbl> 3, 5, 6, 0, 3, NA, 5, 4, 5, 0, NA, 9, 10, 5, 4,...
## $ Enjoy_Wk4      <dbl> NA, 0, 0, NA, NA, NA, 0, NA, 4, 2, NA, NA, NA, ...
## $ Enjoy_Wk8      <dbl> NA, 5, 10, 1, NA, NA, 0, 0, 0, 0, NA, NA, NA, N...
## $ Enjoy_Wk12     <dbl> NA, 4, 10, 5, 10, NA, NA, NA, 4, 0, NA, NA, NA,...
## $ Enjoy_Wk24     <dbl> NA, 4, NA, 4, NA, NA, NA, 10, 5, 4, NA, NA, NA,...
## $ Enjoy_Wk48     <dbl> NA, 8, NA, 7, NA, NA, NA, 5, NA, 4, NA, NA, NA,...
head(bpi)
tail(bpi)

Pain severity scale (PSS)

Clean data

# Generate data
bpi_pss <- bpi %>%
    # Select columns
    select(ID, Site, Group, Sex, 
           starts_with('Worst'), starts_with('Least'),
           starts_with('Avg'), starts_with('Now')) %>%
    # Convert to long format
    tidyr::gather(key = pain_question,
                  value = pain_intensity,
                  -ID, - Site, - Group, -Sex) %>%
    # Seperate pain_question into constituent parts
    tidyr::separate(col = pain_question, 
                    into = c('pain_question', 'time_point')) %>%
    # Calculate PSS
    group_by(ID, time_point) %>%
    mutate(PSS = round(mean(pain_intensity, na.rm = TRUE), 1),
           PSS = ifelse(is.nan(PSS),
                        yes = NA,
                        no = PSS)) %>%
    # Order the time points
    ungroup() %>%
    mutate(time_point = factor(time_point,
                               levels = c('BL', 'Wk4', 'Wk8', 
                                          'Wk12', 'Wk24', 'Wk48'),
                               ordered = TRUE))

Plots

# Create spaghetti plot using all PSS data, to allow for faceting.
# Change time point to numerical so as to space appropriately in graph
# First temporarily convert values from factor to character
bpi_pss <- mutate(bpi_pss,
                  time_point = as.character(time_point))

# Replace values
bpi_pss[bpi_pss=="BL"]<- "BL0"

# Make a copy of bpi_pss 
bpi_pss2 <- bpi_pss

# Separate after 2 characters
bpi_pss <- tidyr::separate(data = bpi_pss,
                           col = time_point,
                           into = c('Text', 'Week'), 2)

# Convert back into numeric form
bpi_pss <- mutate (bpi_pss,
                   Week = as.numeric(Week))

Spaghetti plots

# Create spaghetti plot
# Try using geom.plot, via tutorial at
# https://stats.idre.ucla.edu/r/faq/how-can-i-visualize-longitudinal-data-in-ggplot2/

# Colour by site
ggplot(data = bpi_pss) + 
    aes(x = Week, 
        y = PSS, 
        group = ID, 
        colour = Site) + 
    geom_line()+
    scale_x_continuous(breaks = c(0,4,8,12,24,48),
                       limits = c(0, 48)) +
    scale_y_continuous(breaks = c(0,2,4,6,8,10))

# Colour by individual, facet by site
ggplot(data = bpi_pss) +
    aes(x = Week,
        y = PSS,
        colour = ID) +
    geom_line() +
    facet_grid(Site ~ .) +
    labs(title = 'Pain Severity Scale', 
         subtitle = 'Faceted by study site, coloured by ID')+ 
    theme(legend.position = 'none') +
    scale_x_continuous(breaks = c(0,4,8,12,24,48),
                       limits = c(0, 48)) +
    scale_y_continuous(breaks = c(0,2,4,6,8,10))

# Facet by study site and gender, colour by ID
ggplot(data = bpi_pss) +
    aes(x = Week,
        y = PSS,
        colour = ID) +
    geom_line() +
    facet_grid(Site ~ Sex) +
    labs(title = 'Pain Severity Scale', 
         subtitle = 'Faceted by study site, coloured by ID')+ 
    theme(legend.position = 'none')+
    scale_x_continuous(breaks = c(0,4,8,12,24,48),
                       limits = c(0, 48)) +
    scale_y_continuous(breaks = c(0,2,4,6,8,10))

# Facet by ID and U1 only
bpi_pss_U1 <- bpi_pss %>%
    filter(Site == "u1")

ggplot(data = bpi_pss_U1) +
    aes(x = Week,
        y = PSS,
        colour = ID) +
    geom_line() +
    facet_wrap(facet = ~ID) +
    labs(title = 'Pain Severity Scale', 
         subtitle = 'U1 site only; faceted by ID')+ 
    theme(legend.position = 'none')+
    scale_x_continuous(breaks = c(0,4,8,12,24,48),
                       limits = c(0, 48)) +
    scale_y_continuous(breaks = c(0,2,4,6,8,10))

# Facet by individual at each site
## U1
pss_spagID <- bpi_pss %>%
    select(-pain_question, -pain_intensity) %>%
    unique(.) %>%
    group_by(Site, Group) %>%
    nest() %>%
    mutate(Number = row_number()) %>%
    unnest() %>%
    nest(-Number) %>%
    mutate(plot = map(data, 
                      ~ggplot(data = .x) +
                          aes(x = Week,
                              y = PSS,
                              colour = Sex) +
                          geom_line() +
                          geom_point(size = 2) +
                          labs(title = 'Pain severity scale',
                               subtitle = glue('Site: {str_to_title(unique(.x$Site))} ; Group: {.x$Group}')) +
                          scale_x_continuous(breaks = c(0,4,8,12,24,48),
                                             limits = c(0, 48)) +
                          scale_y_continuous(breaks = c(0, 2, 4, 6, 8, 10),
                                             limits = c(0, 10)) +
                          facet_wrap(facets = ~ID, ncol = 4) +
                          theme(legend.position = 'top')))

pss_spagID$plot
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

Summary plots

## PSS: Faceted by study site, coloured by group
ggplot(data = bpi_pss2) +
    aes(x = time_point,
        y = PSS,
        colour = Group,
        fill = Group) %>%
    geom_boxplot(alpha = 0.6) +
    facet_grid(Site ~ .) +
    labs(title = 'Pain Severity Scale (PSS)', 
         subtitle = 'Faceted by study site, coloured by group')

## PSS: Faceted by study group, coloured by sex'
ggplot(data = bpi_pss2) +
    aes(x = time_point,
        y = PSS,
        colour = Sex,
        fill = Sex) %>%
    geom_boxplot(alpha = 0.6) +
    facet_grid(Group ~ .) +
    labs(title = 'Pain Severity Scale (PSS)', 
         subtitle = 'Faceted by study group, coloured by sex')

## PSS: Faceted by study site and group, coloured by sex
ggplot(data = bpi_pss2) +
    aes(x = time_point,
        y = PSS,
        colour = Sex,
        fill = Sex) %>%
    geom_boxplot(alpha = 0.6) +
    facet_grid(Site ~ Group) +
    labs(title = 'Pain Severity Scale (PSS)', 
         subtitle = 'Faceted by study site and group, coloured by sex')

## PSS: Women only
ggplot(data = filter(bpi_pss2, Sex == 'female')) +
    aes(x = time_point,
        y = PSS) %>%
    geom_boxplot(alpha = 0.6) +
    labs(title = 'Pain Severity Scale (PSS)',
         subtitle = 'Women only')

## PSS: Women only (excluding U1)
ggplot(data = filter(bpi_pss2, Sex == 'female' & Site != 'u1')) +
    aes(x = time_point,
        y = PSS) %>%
    geom_boxplot(alpha = 0.6) +
    labs(title = 'Pain Severity Scale (PSS)', 
         subtitle = 'Women only (excluding u1)')

Pain interference scale (PIS)

Clean data

# Generate data
bpi_pis <- bpi %>%
    # Select columns
    select(ID, Site, Group, Sex,
           starts_with('Activ'), starts_with('Mood'),
           starts_with('Walking'), starts_with('Work'),
           starts_with('Rel'), starts_with('Sleep'),
           starts_with('Enjoy')) %>%
    select(-contains('Relief')) %>%
    # Convert to long format
    tidyr::gather(key = interference_question,
                  value = interference_intensity,
                  -ID, -Site, -Group, -Sex) %>%
    # Separate pain_question into constituent parts
    tidyr::separate(col = interference_question, 
                    into = c('interference_question', 'time_point')) %>%
    # Calculate PIS
    group_by(ID, time_point) %>%
    mutate(PIS = round(mean(interference_intensity, na.rm = TRUE), 1),
           PIS = ifelse(is.nan(PIS),
                        yes = NA,
                        no = PIS)) %>%
    # Order the time points
    ungroup() %>%
    mutate(time_point = factor(time_point,
                               levels = c('BL', 'Wk4', 'Wk8', 
                                          'Wk12', 'Wk24', 'Wk48'),
                               ordered = TRUE))

Plots

# First temporarily convert values from factor to character
bpi_pis <- mutate (bpi_pis,
                   time_point = as.character(time_point))

# Replace values
bpi_pis[bpi_pis=="BL"]<- "BL0"

# Make a copy of bpi_pis
bpi_pis2 <- bpi_pis

# Separate after 2 characters
bpi_pis <- tidyr::separate(data = bpi_pis,
                           col = time_point,
                           into = c('Text', 'Week'), 2)

# Convert back into numeric form
bpi_pis <- mutate (bpi_pis,
                   Week = as.numeric(Week))

Spaghetti plots

# Colour by site
ggplot(data = bpi_pis, aes(x = Week, y = PIS, group = ID, colour = Site))+ 
    geom_line()+
    scale_x_continuous(breaks = c(0,4,8,12,24,48),
                       limits = c(0, 48)) +
    scale_y_continuous(breaks = c(0,2,4,6,8,10))

# Colour by individual, facet by site
ggplot(data = bpi_pis) +
    aes(x = Week,
        y = PIS,
        colour = ID) %>%
    geom_line() +
    facet_grid(Site ~ .) +
    labs(title = 'Pain Interference Scale', 
         subtitle = 'Faceted by study site, coloured by ID')+ 
    theme(legend.position = 'none') +
    scale_x_continuous(breaks = c(0,4,8,12,24,48),
                       limits = c(0, 48)) +
    scale_y_continuous(breaks = c(0,2,4,6,8,10))

# Facet by study site and gender, colour by ID
ggplot(data = bpi_pis) +
    aes(x = Week,
        y = PIS,
        colour = ID) %>%
    geom_line() +
    facet_grid(Site ~ Sex) +
    labs(title = 'Pain Interference Scale', 
         subtitle = 'Faceted by study site, coloured by ID')+ 
    theme(legend.position = 'none')+
    scale_x_continuous(breaks = c(0,4,8,12,24,48),
                       limits = c(0, 48)) +
    scale_y_continuous(breaks = c(0,2,4,6,8,10))

Summary plots

## PIS: Faceted by study site, coloured by group
ggplot(data = bpi_pis2) +
    aes(x = time_point,
        y = PIS,
        colour = Group,
        fill = Group) %>%
    geom_boxplot(alpha = 0.6) +
    facet_grid(Site ~ .) +
    labs(title = 'Pain Interference Scale (PIS)',
         subtitle = 'Faceted by study site, coloured by group')

## PIS: Faceted by study group, coloured by sex'
ggplot(data = bpi_pis2) +
    aes(x = time_point,
        y = PIS,
        colour = Sex,
        fill = Sex) %>%
    geom_boxplot(alpha = 0.6) +
    facet_grid(Group ~ .) +
    labs(title = 'Pain Interference Scale (PIS)',
         subtitle = 'Faceted by study group, coloured by sex')

## PIS: Faceted by study site and group, coloured by sex
ggplot(data = bpi_pis2) +
    aes(x = time_point,
        y = PIS,
        colour = Sex,
        fill = Sex) %>%
    geom_boxplot(alpha = 0.6) +
    facet_grid(Site ~ Group) +
    labs(title = 'Pain Interference Scale (PIS)', 
         subtitle = 'Faceted by study site and group, coloured by sex')

## PIS: Women only
ggplot(data = filter(bpi_pis2, Sex == 'female')) +
    aes(x = time_point,
        y = PIS) %>%
    geom_boxplot(alpha = 0.6) +
    labs(title = 'Pain Interference Scale (PIS)', 
         subtitle = 'Women only')

## PIS: Women only (excluding U1)
ggplot(data = filter(bpi_pis2, Sex == 'female' & Site != 'u1')) +
    aes(x = time_point,
        y = PIS) %>%
    geom_boxplot(alpha = 0.6) +
    labs(title = 'Pain Interference Scale (PIS): Women only',
         subtitle = 'Excluding U1')

Session information

sessionInfo()
## R version 3.4.2 (2017-09-28)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Sierra 10.12.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] bindrcpp_0.2  ggplot2_2.2.1 glue_1.1.1    stringr_1.2.0 tidyr_0.7.2  
## [6] purrr_0.2.4   dplyr_0.7.4  
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.13     knitr_1.17       bindr_0.1        magrittr_1.5    
##  [5] hms_0.3          tidyselect_0.2.2 munsell_0.4.3    colorspace_1.3-2
##  [9] R6_2.2.2         rlang_0.1.2      plyr_1.8.4       tools_3.4.2     
## [13] grid_3.4.2       gtable_0.2.0     htmltools_0.3.6  lazyeval_0.2.0  
## [17] yaml_2.1.14      assertthat_0.2.0 rprojroot_1.2    digest_0.6.12   
## [21] tibble_1.3.4     reshape2_1.4.2   readr_1.1.1      evaluate_0.10.1 
## [25] rmarkdown_1.6    labeling_0.3     stringi_1.1.5    compiler_3.4.2  
## [29] scales_0.5.0     backports_1.1.1  jsonlite_1.5     pkgconfig_2.0.1